FINAL REPORT 


{jyo A-* /t/2 


///-O' /' ^ 

=25 3 


3if 


TO 

THE NATIONAL AERONAUTICS AND SPACE AGENCY 

FOR 

POWER SPECTRAL ESTIMATION ALGORTHMS 
(CONTRACT No. NAG-5-499) 


MANJIT S. BHATIA 
PRINCIPAL INVESTIGATOR 

BOWIE STATE UNIVERSITY 
BOWIE, MARYLAND 

DECEMBER, 1989 


(, ASA-CR-ls^jj 
ALbORIfn^^ p |r 
34 o 


POWER 

Report 


SPEC TRAL 
( bow i e 5 


fstim ation 
t '*te Col 1 . ) 
c<l CL 09B 


0->/ 61 


*90-1401^ 

Unc las 
02b3Q2> 


\ 

\ 

\ 


ORIGINAL PAGE IS 
OF POOR QUALITY 


POWER SPECTRAL ESTIMATION ALGORITHMS 


In — this project, we have developed 
algorithms to estimate the power spectrum. . 
using Maximum Entropy Methods;”' These 
algorithms were coded in FORTRAN 77 and were 
implemented on the VAX 780. 

We recall that the important 
considerations in this analysis are: 


,'l) RESOLUTION:',,, How close in frequency can"' 
two spectral components be spaced and still 
be identified?^- 

*2^ .DYNAMIC RAN&El t/r Hbw small carT'a spectral 
peak x be, relative to the largest, and still 
be observed in the spectra? 

3> VARIANCE:. How accurate is~ the estimate of 
the spectra to the actual spectra? 

^Dur work with the application of the 
algorithms based on Maximum Entropy Methods 
to a variety of data shows that these 
criteria are met quite well. Additional work 
in this direction would help confirm our 7 
findings . 

All of the software developed has been 
turned over to the technical monitor. A copy 
of a typical program is included in™ this 
report . 

Some of the actual data and graphs based 
on this data are also included in ^this 
report . 



2—D Maximum Entropy Method 

CASE 1 : DFT length = 8 



(Xj'xj)s 


Index 



> 

v CASK. \ » 


NUMBER UH SINUSOID?) - 1 
NOISE- POWER « 5.0 


SINUSOID = 1 

POWER = 1. . 0 
XFREQ - 0.375 
YFREQ - 0.. 25 

ORIGINAL PAGE IS 
OF POOR QUALITY 
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2—D Maximum Entropy Method 

CASE 2 : DFT length = 16 
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2—D Maximum Entropy Method 

CASE 3 : DFT length — 32 
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0 n 03 1 30 

0 „ 69 

■ h 6 

450 

0 .. 43750 

0 , 06250 

0 7 E 

i. 

45 1 

0 - 43750 

0 n 09380 

0 „ 9 :... 

» / / 

452 

On 43750 

0 „ i 2500 

1 1 i:i 

>'6 

453 

On 43750 

0. .1 5630 

1 n 32 

i / 

454 

0 > 43750 

On 18750 

2 „ 20 

*20 1 

455 

0 n 43750 

0 2 1 880 

y Q 

: J : %: ■“> 

456 

0 * 43750 

O „ 25000 

. 

. .. » A / 

457 

U « 43 / 50 

0 n 28 1 30 

o r:/ 
„ y % ,.f 

524 

458 

0 „ 43750 

0 .. 3 1 250 

2 : .. 20 

20 1 

459 

0 n 43750 

0 n 34380 

1 n 58 

72:1 

460 

0 . 43750 

0 „ 37500 

] -t . *i n; ,-j / 

.1 . i h»J 70 

461 

0 « 4 3 7 5 0 

0„ 40630 

On 93 

/ / 

462 

On 43750 

0 , 43750 

On 78 

6 1 5 

463 

On 43750 

O ,, 46880 

0 69 

/;;> . 3 

464 

0 . 43750 

0 n 50000 

0 n 64 

878 

465 

On 43750 

0 - 53 1 30 

On 63 

286 

466 

0 43750 

0 „ 56250 

On 64 

- J, 

467 

On 43750 

0 .. 59380 

0 n 67 

717 

468 

0 . 43750 

0 . 62500 

0 „ 73: 

2 5 6 

469 

0„ 43750 

0 n 65630 

On so: 

377 

470 

. 0 „ 43750 

0 . 68750 

On 87’ 

990 

47 ;j 

0 n 43750 

0 „ 7 1 880 

On 94* 

094 

472 

0- 43750 

On /bo 00 

0 n 96- 

4 6 9 

473 

0 „ 43750 

0 n 78 1 30 

O , f 9 '1 *' 

794 

474 

0.. 43750 

0 . 8 1 250 

0.. 87' 

790 

475 

On 43750 

0 - 84380 

On so; 

'7 ' r 

476 

0 n 43750 
On 43750 

r .) „ c> 7 5 J y ) 

On 73, 

235 

4 7 7 

0 „ 9067.0 

O .. 6 / 

1 j 
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' 47 B 

0 . 43750 

0. 93/ bC. 


479 

0. 43750 

0 .. 9688 V 

C ) „ r\ , ■ > 

480 

0 „ 46880 

0 „ 0000*0 

t J r, r. } " r ■■ !■ • /I - * 

48 .! 

On 46880 

On 031 30 

O. 6 1 \ 

482 

0 „ 46880 

0 .. 06250 

On 34360 

483 

0 n 46oo'0 

0 tl 0 9380 

0 « 81 ; ’ i 

484 

■ 0 „ 4 68 SO 

Or 12500 

1 r 03827 

485 

0.46880 

0. 15630 

;! 3013 !!. 

486 

0. 46880 

Or, 18750 

1 659 1 7 

487 

0 ,, 46880 

0 » 2 1 880 

2 „ 02389 

488 

0. 46880 

0. 25000 

2, 19245 

4S9 

0. 46880 

0 * 28 1 30 . 

2r 02389 

490 

0. 46880 

0 >■ 3 1 250 

1 * 659 1. 7 

49 1 

0. 46880 

0., 34380 

1 . 3043 :l 

492 

0 . 46880 

u „ 3 /boO 

.1 „ 03827 

493 

0. 46880 

0„ 40630 

0 „ 85840 

494 

(j „ 4 6^ 8 8 O 

0. 43750 

0 „ 74360 

495 

0.46880 

0.46880 

(..) „ 6 / 6Ub 

496 

0. 46880 

On bOUUO 

0 » 64439 

497 

0 . 46880 

0 53 1 30 

0 „ 64261 

498 

0« 46880 

0 rr 56250 

U •> <9 6 6 b 6 

499 

Or, 46800 

0 . 59380 

Or 72299 

500 

0. 46880 

0 „ 62500 

0 ,, 803 1. 2 

501 

0 „ 46880 

0 .1 65630 

„ v :! 353 

502 

0. 46880 

0. 63750 

1. « 03399 

503 

0.46880 

0 - 7 1 880 

1 « 1 35 1 2 

504 

0 » 46880 

0- 75000 

l .. 1 / 5/ j 

505 

0 .46880 

0 „ 78 1 30 

1 „ 13512 

bo 6 

CO 

CO 

-0 

's 

0 . 8 1 250 

1 „ 03 390 

507 

0. 46880 

0- 84380 

0 i( 9 ;| 353 

508 

0. 46080 

Or 87500 

On 805 :i i 

509 

0 . 46880 

0 90630 

0.. 71 '24 9 

5 1 0 

0. 46880 

0 .. 93750 

0 n 6' 68 56 

5 :!. l 

0» 46880 

0 .. 96880 

(..) „ 64261, 

512 

0 . 50000 

0,. 00000 

On 64292 

5 :L 3 

0 „ 50000 

0 03 1 30 

0 n 6575 1 

5 1 4 

0 * 50000 

O „ 06 2' 50 

0 70293 

5 1 5 

0, 50000 

0.. 09380 

0„ 78398 

5 .1 6 

0 „ 50000 

Or. 12500 

0 . 90 740 

517 

0. 50000 

0., 15630 

1 ■ 0 7693 

5 1 8 

0. 50000 

On 18750 

1 27961 

5 1 9 

0. 50000 

0 - 2 1 880 

1 n 46366 

520 

0. 50000 

0 „ 25000 

1 „ 54 1 52 

52 1 

0» 50000 

0 rr 28 1 30 

1 n 46366 

'""i <r-l 

si. .e. 

0. 50000 

0 „ 3 1 250 

1 „ 27961 

523 

0 „ 50000 

0 „ 34380 

1 ,t 07693 

524 

0 « 50000 

0 * 37500 

0 „ 90740 

525 

0 50000 

0.. 40630 

0,. 7 83 Vo 

526 

0. 50000 

0 M 43750 

0„ 70293 

527 

Or 50000 

On 46880 

j.--, / c:r Trr -t 

V « £> , J / ,.J t 

528 

0 « 50000 

0 „ 50000 

o „ ooooo 
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2—D Maximum Entropy Method 

CASE 4 : DFT length = 10 





Index 


• CASE 4: 

NUMBER OF SINUSOIDS 
NOISE POWER 6.0 


SINUSOID 


1 

POWER 

= 

1 .. 0 

XFREQ 

= 

0 .. 1 

YFREQ 

zz 

0 , . 1 . 

SINUSOID 


.fL 

POWER 

r=: 

1 - 0 

XFREQ 


0 . 3 

YFREQ 

= 

0- ;i 

SINUSOID 

ssr 

3 

POWER 


1 0 

XFREQ 

= 

0 2 

YFREQ 

s= 

0 „ 2 


ACF MATRIX - 7 >: 


DFT LENGTH = 10 
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INDEX 

Fm 

Fy 

s <F;, , Fy 

0 

Om ouuoo 

0 „ 00000 

0. OOOOO 

1 

0 n 00000 

Or 10000 

1 iv 05745 

T «* 

0 . 00000 

0 n 20000 

0 . 6 A 67 1 


0 „ 00000 

0 » 30000 

0 n 50208 

4 

0 . 00000 

0- 40000 

0 . 6 1 538 

\\z 

0 „ 00000 

0„ 5 GOOD 

0. 8:1. 94 o 

6 

0 n 00000 

0 60000 

0 . 00000 

7 

0 .. 00000 

0 if 70000 

0 . OOoOO 

S 

0 „ 00000 

o . auoou 

0 . OOOOO 

9 

0 r 00000 

0 .. 90000 

0 n OOOOO 

10 

0. 10000 

0 « 00000 

0 r 73649 

1 1 

0 . 1 0000 

0. 10000 

6 .. 40696 

1 2 

0 » 1 0000 

0 .20000 

1 . 13217 

1 3 

0.. 10000 

O m 30000 

0„ 59657 

14 

0. 10000 

0 » 40000 

0 ,i 73546 

15 

0 10000 

0 « 50000 

0„ 66967 

16 

0„ 10000 

0 w 60000 

0.. 5 BO 1 2 

17 

Or 10000 

0 r 70000 

o r a 1 9 :?:; i 

18 

Or 10000 

0 . 80000 

0 R 65055 

19 

Or 10000 

o „ 90000 

Or 43186 

20 

0 if 20000 

0 ii 00000 

0 .. 6496 / 

21 

0 r 20000 

o, i on on 

2« 13862 

v:. s... 

0 » 20000 

0 rv 20000 

6 r 19231 

..... 

Or 20000 

0 h 30000 

1 „ 04869 

24 

0 . 20000 

0„ 40000 

Or 62238 

25 

0- 20000 

0„ 50000 

0 .• 6 1 939 

26 

Or 20000 

0 „ 60000 

0 R 696 1 3 

..... -j 

0 n 20000 

Or 70000 

0 / 4 5 6 4 



28 

0 «i 20000 

0 , 80000 


... '' 1 


• • 

29 

0 • 20000 

0„ 90000 

i' > K 

2 


30 

30 

0 „ 30000 

0 „ 00000 

1 . 

i 


-■ f 

31 

0 „ 30000 

0.. 10000 

4 .. 

■ 



S 2 

0 * 30000 

D . OOi H K ) 

1 „ 

O 

;| 

03 


0 „ 30000 

0 n 30000 

0 

37 


0 .6 

34 

0 .. 30000 

0 .. 40000 

< ) u 

?■ , i 

i::' 

•'i ■ ; 

■! , t 

35 

Oh 30000 

0 ir 50000 

0 H 

7 7 


7' / 

36 

0 ■ 30000 

0 , 60000 

0 H 

6 1 

3 

Vz< 

37 

0„ 30000 

0„ 70000 

0 

r.; ' i 

».,* / 


33 

38 

0 „ 30000 

0 » 80000 

0 , 

/ 

8 

3 :i 

3 v 

0 30000 

0 „ 90000 

a „ 

/ , j 

7 

r. . ■ > 

40 

0 „ 40000 

0 . 00000 

0 . 

69 

3 

0-1 

4 1 

0„ 40000 

Oh 10000 

.1- V 

60 

r. , 
-..J 

49 

42 

0„ 40000 

0 w 20000 

0 H 

56 

Q 

48 

43 

0 „ 40000 

0 H 50000 

0 h 

43 

9 

!::) 5 

44 

0 , 40000 

0 - 40000 

0 .. 

77 


85 

43 

0 „ 40000 

0 .. 50000 

0 H 

/ ".i 

i 

. 1 . 

43 

46 

0 « 40000 

0. 60000 

On 

in; 

4 

r ) / 

47 

0 „ 40000 

On 70000 

On 

83 

2 

45 

48 

0 „ 40000 

0 » 80000 

Kj „ 

75 

0 

:! 2 

49 

0 - 40000 

0 » 90000 

I } 

43 

o 

2 4 

50 

0„ 50000 

0 « 00000 

0 H 

43 

y 

03 

51 

0 n 50000 

o. :i oooo 

() n 

46 

8 

i..,) 

t" *■-. 
j x* 

0 . 5 DC) 00 

0 „ 20000 

f") 

SB 

'{ 

. 1 . 6 

3 5 

0 50000 

0 .. 30000 

\) 

/ ' y 

■i 

69 

54 

0 „ 50000 

0 h 40000 

0 H 

67 


7 4 

rs* r~ 
JJ 

0 „ 50000 

0 .. 50000 

0 U 

C) j) j 


Q ( ) 
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RESULT Witll NEM 




PROGRAM MR II 2D 


c 'K‘ ^ X- *>f ‘K‘ X- H* A- ‘ff -X- K- -K- H- -)f -X- K* If v -X X- if Yc -X- if if -rr -if if *K -if *K* -rf >:■ if «• If -fr i •■ Of X •£ H- K- if *K Of if if if ¥. X- ->f K- if if if -K^ *X- if If X- -If if 


c 

c 

c 

c 

c 

c 

c 

c: 


THE 2-D MAXIHUM ENTROPY METHOD PSD ESTIMATOR: 

1 o ob t a ;L n t he 2—D m e m s a I li l; i o n f a r r e a I 
d at a ( r e a 1 s y rr \ rn e t r i c a c: f ) , u s i n g d :i. r e c t cl t l: 
computation. The MEM PSD does not required that 
the k n o w n 1 a g s o f t hi e 2 — D a c f h a v e a li n :i f o r m s u p p o r t 

gi id? arbitrary lags can be used as the constraints. 


c 

x I am 

=s 

I a m b cJ a c a e f f i c i e n L a r r a y . 

c: 

x o 1 d 


old la in b c j a a r r a y requir e cl -f : o r 

c 



beta computed: i an. 

c 

r 

= 

k n o w n a c t v a 1 u e . 

c 

px 


computed act values or correction acf . 

c 

gap 

rr 

logic a 1 a r r a y s p e c i y i n g w h ei t g a p s 

c 



i f any, exist in the* ac*f . 

c 

n 

= 

df t 1 eng th used i n the i ter at i ons. 

c 

r,2 


n/2 

c 

n 2 1 

= 

n 2 + 1 

c 

n4 

- 

n -X'- n 

c 

n 1 

= 

in ax + ve ( x , y ) i n ri ex f or i: h e smal 1 est 

c 



s <.j u <_-< r (..• r e g i o n c o ri t a .i n i n g *' h; ' „ 

c 

inn 


m a >: a c -f* a r r a y s i :: e for t h i s p r o g r a rn ( 2 5 

c 

mn2 


center point of the acf array. 

c 

mn3 

rr; 

mn2 + 1 

c 

mn4 

-- 

2 * mn2 

c 

mi 

= 

rnin (x,y) index far- known region. 

c 

m2 

rr: 

max (Xj,y) index for known region. 

c 

sc:l t 

rr: 

scale f a c t o r . 

c 

:l:st 

= 

error 1 evel . 


) . 


c *K* * * if if if if * Of if if Of -X- X- if Of if if * if -X- K- X- Of Of Of Of Of if Of Of -X- Of Of -X X- Of Of ¥. K X- Of if -If Of if Of Of -X- -If Of ->f -X* ■ 


real *4 xl, 

\* S' v 

5 j - «i 

den 





real *4 x 1 am 

(25,25) 

, xol d ( 

»:r 

.ClwJ, 

O nr \ tri' 

U ) f f r . , 

25), px (25 

common n , n2 

, n 2 1 ,, 

n 4 , 

nl. 

mn 

? m n 2 , ni n 3 , m n 4 , m 1 , 

format ( lOx , 7 

n i t r 

? 

i 3, 

5x , 

? error “ 

7 , e J 2 . 4 , / ) 

format ( lOx , 7 

a 1 p h a 

? 

e 1 4 


/) 


format < i Ox , 7 

beta 

’ ;t 

e 1 4 

. 6, 

/) 


mn = 25 







mn2 =5 (mn - 

1 ) / 2 

+ i 





call acf 2d (r 

) 
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/ rum 2 


rmn 



r ( 

mn2 ? mn2 ) 

cJ o 

i -- 

m 1 

m2 

do 

j ~ 

ml 

? m2 


r(i, 

j > 

- r (i ? .j ) 


e n ddo 


. « > r- \ m 


end do 


old = 1 . 0E30 
zoldl = 1.0E30 


write!*, 90) 

? 0 far nr at ( / , I 0 x , 7 e n t e r e r r a r 1 a v e 1 ( r: t *:> i ) 7 

read 2 tst 

wri te <■*, 95 ) 2 tst 

?5 formate/, 5x , 7 error level ( 2 tst) = 7 , f 

self — 0.5 
alpha = 0.0 
beta - 0.0 

do i = 1 , nin 
do j r - 1 , inn 

x o 1 d ( i , j ) = 0 - 0 
x 1 am ( i , j ) - 0. 0 

p >: C i , j ) - 0. 0 

en ddo 
enddo 

x 1 am (mn2, mn2) = ;[ . 0 

>:ol d (mn2, mn2) - xl am (,mn2, mn2) 

den = 0.0 
do i — ml, m 2 
do j - ml, m2 

den -- den + r ( i , j ) v r ( i , j ) 
enddo 
en ddo 

den = den - 1.0 
n i t r ~ 0 

>00 cal 1 ft Cpx , xl am, x 1 , x2) 

nitr « nitr + 1 
z ~ 0 . 0 
do i = ml, m2 

do j =s ml, m2 

px (i , j j = r (i , ;j ) - px (i , j > 

z - z -1- p x ( i , j > ■* p x ( i , j ) 

enddo 
end do 


Z •— 2 

/ den 

wr i t e 

( * , 3 . ) nitr, z 

if (2 

•It. 2 tst ) go to 9000 


/ ) 


U . 8 , / ) 
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if 


( alpha . gt. O - 9 V 9 9 ) go to 9000 
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i r ( z .ge. z c 1 d .or. b ot a . n e . 0 .. 0 ) goto 970 
go to 971 


970 alpha = ( 1.0 + alpha ) / 2.0 
sc It = self / 2.0 

z c j 1 d 1 = z o 1 d 

971 zold = z 


c a 1 J f 1 2 ( p x , >: 3 ) 



>; X ~ 

1 . 0 

+ self 

* x 2 

/ x 3 


a 1 p h a 

= a max 1 ( 

a ]. p h a 

, XX ) 


wr i t e 

( * 

<1 +L ) ct 

1 pha 



cal 1 

-ft 3 

( ;•! o 1 d 

V P X * 

a 1 pha, x 1 a 


call 

f t2( 

>: 1 am. 

x 4 ) 



if ( : 

x 4 . 

gt. 0. 

0 ) g o 

to 972 


b m i n 

=: — 

: : 4 / 

( >! 1 - 

x 4 ) 


beta 

= ( 

1 . 0 + 

( 1 . 0 

- self ) 


al phi a 

( 

1 . 0 + 

a J. pha 

) / 2 . 0 


g cd t o 

973 




?72 

beta = 

- 0. 1 

0 



?7 3 

wr i te 

( * 

n 3 ) beta 



d o i : - 

= ml 

m2 




d o j = 

= ml j 

i m2 




xlam( i , j ) — beta * MOLD ( i , j ) 
Mold (i , j ) ~ x 1 am ( i , j ) 

enddo 
en d do 

go to 1 000 
*000 conti nue 


1.0 / b iTi i n — i . C> 


+ f 1.0 — beta ) -ft 


write(^ ? 9900 ) 

9 U O format ( lx, r con v e r g e n c e a t: hi :i e v e? d ! t e s t i n p r og r e s s 
* pat i ence. ... ? ) 


n 

= 4 

■a- 

n 

n 2 

” n 

/ 


n 21 

= n 2 

+ 

1 

n 4 

= n 

K- 

n 

c: a 1 ]. 

ft < 

p;: 

1 5 


) * Li m i n 


x 1 am ( i , j ) 
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■z ~ 0 . 0 
do i = ml, m2 
do j =■- ml, m2 

- ~ 2 + ( r ( i , j > - px ( i , j ) > **• 2 

enddo 
enddo 

z = z / den 

if ( z . 1 e. ztst ) go to 974 



write ( *, 

9002 ) 

z 


9002 

format < / , 

? error 

■ level too h 

i i g hi : 


go to 975 




974 

wr i te ( x. 

9001 ) 

z 


9001 

f or mat ( / , 7 

good E 

solution, err 

or i s 

975 

cent 1 nue 




cc 

do i ml 

? m2 



cc 

do j -= ml 

5 m2 



cc 

write < x. 

9300 ) 

i ? j .| x 1 am ( 

i ? J 

9300 

format ( :1 k 

, ’ lambda ( ’ , 12, IX., 

13, 5 ) 

cc 

wri te f 7, 

9400 ) 

x 1 am ( i , j 

) 

9400 

f g i" m a t ( 1 x , 

f 10.7 

) 


cc 

enddo 




cc 

enddo 





c a 11 pr i n t ( 

x 1 am ,, 

m2-ml+l ) 



e 12 . 4 ,// ) 


, E14.6, /) 


stop 

end 


G * -K * * X X X X * Yr X X -X- >; * K X X X X X X X X K X -X- * -X- * * * * X- X X * X X if -X- X -‘X- X X -X- -X X -X- -X- X X- -X* X -X- X X -if X X X X X X 

sub r out i n e ac f 2d ( r ) 

(2 *X* -X- X X' X X X -X- X -X -X- X -X* -X* K* -X X* X* X* -X* X X X *}f -if X X -if -X X * X -X- X X X X X X X X X X X- X X -X X -if X X- X X X X X X X X X X X X X X 


i n teg or 

n ., n 1 , run i , 

mr»2, 

iTl :l 

m2 


real *4 

r (25, 25) , 

p (20) 

? *'* 

freq ( 20) , 

y f req ( 20 ) , noi se 

common 

n, n 2 , n21 

•I n 4 , 

n 1 

mn, mn2ji 

mn 3 9 mn4, ml, m2 


pi ~ 4.0 * at an (1.0) 
wri te (x, 110) 

1 10 format </, lOx , ? enter number of sinusoid s ? , /) 

read*, nsi n 
w r i t e ( * , l ;(. 2 ) n s i n 

112 format(/ t 5:-:, ? number of sinusoids ~ ", i 2, /) 



do i 


1, ns in 
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130 

wr i te ( * , loO) 
•format ( / , 1 Ox 


read * ? p ( i ) ? 

132 

write ( 132 

format ( / , 5x , 
* ? > ; f r e q = 7 ,, 


en ddo 

150 

write ( *, 150) 

f or mat ( / < 1 Ox 


read-*., noine 

152 

write( * ? 1 52 
format ( / 5x ? 

170 

wr i. t e ( * 7 1 70 ) 

format ( / ? lOx 


power = 


f 7 . 3 , 


f 8. 


* acf (must be odd)’,/) 


172 


19: 


199 


read^-j 

n 12 


write ( 

K-, 172 

) n 1 2 j, 

f or mat ( / , 5x , 

? act 

11 

r— 1 

c 

( nl2 - 

1 ) / 

ml = 

mn2 — n 1 


rri2 = 

mn2 * nl 


mn3 = 

mn2 + 1 


mn4 = 

2 *■ mn 


wr i t e ( * , 19/) 

format (/ , 10::, 

7 df t 

r e a ci * , 

n 


wr i te ( 

•*, 199 

) n 

format ( / ? 5x , 

7 df t 

r.2 « 

n / 2 


n21 = 

n2 + 1 


n4 = 

n * n 


do 30 

ns = 1 , 

nsi n 


wx 2 - 

*■ pi 


wy = 2. 

•* p i 

do 30 

i = rn 1 1 

m2 


ia = i - 

■ on n 2 

do 30 

j = m 1 ? 

m2 


ja = j - 

- mn2 


, i 2 , ’ BY i 2 , / ) 



if ( ns . eq. 1 ) r < i , j ) =-0.0 

r ( i , j ) = r ( i , j ) (ns) -k-cos (wx**i a+wy*j a) 


30 continue 


r( mn2 ? mn2 ) = r ( mn2 ? rnn2 ) + noise 

DO I = * M 1 , M2 

WRITE (*, 91 .1 ) ( R(I,J), J“M 1 , M2 ) 

91 1 FORMAT ( 9 <F:L0. 4) ) 

ENDDO 

return 

end 

C -H-*****^ Jf*^*‘K *^^ )f**^*Jf^^ f{ ^**^**^-^*>: *^ * **^ ^#-)f* K- -K* 4r ¥s ■#■ K* X -ft * -K- # K* -X- 45- * * X- -if 
subrout i no f t2 (px ? x mi n ) 

real *4 p>:(25,25) 

comp ]. ex z0 ? zl ? 22 (1 z 3 fl z 4 3 z 5 ? :■: 6 ? z 7 

common n ? n2 ? n21 , n4 , n I , mn ? m n 2 ,, mn3, mn4 , m 1 , m2 

pin = 8 . 0 x-at n;i ( 1 . O') / T 1 oa t ( n ) 

x 00 = 0.0 
x 22 = 0.0 
do 100 i = ml, m2 
da 100 j = ml, m2 
x 00 = x 00 + p x ( i , j ) 

1 00 x 22 = x 22-i-px ( i , j j * ( (-1 . 0 ) *x < i + j -mn4 ) ) 
x mi n = ami n 1 ( x 00 3 x 22 ) 

x 0 = p x ( mn 2 ? mn 2 ) 
z 0 = c mp 1 x ( c os ( p i n ) , si n ( p i n ) ) 
z 1 = (1.0, 0 . 0 ) 

z6 = cmp 1 x (cos ( pi n*n 1 ) ? -si n (pin *n 1 ) ) 
z 7 = (1.0,0. 0 ) 


200 


do 600 k = l ? n21 
z 2 = z 1 
x 1=0. 0 

do 20 0 m = m n 3 , m 2 

x 1 = x 1+px <rn, mn2) ■* (real (z2) ) 
z2 = z2*zl 


300 


z 3 = ( 1 . 0 ? 0 . 0 ) 

do 300 1 = l H n 

i -F (k*l -eq. 1) go to 500 

if (k. eq. 1 .and. 1 .gt. n21) go to 601 
i -f (k. eq. n21 .and. I .gt. n2) go to 600 
x 2 = 0. 0 
z 4 = z 7 

do 400 m = ml, m2 
z 3 — z o 

do 300 nn = mn3, m2 

>: 2 = x 2 + p x ( m , n n ) x- r e a 1 ( z 4 x- z 5 ) 
z 5 = z 5#z 3 
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1 t ( rn . e q . m2 ) go 1: o 4 00 
2 4 = 24*21 

400 cant i nue 

>; k I - x 0 + 2 • 0 * ( x 1 + x 2 ) 

:>: mi n = ami n 1 ( x mi n , >: k I ) 
500 = 2 3*r.O 

601 zl ~ 2 1*2 0 

2 7 = 2 6 *c on j g ( 2 4 ) 

600 conti nue 
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return 

end 

c * * * * * * ******** * * ■¥: * K * * * * * * * -V:- * -X- * * * * * * * * * * * * * * H X * * * * * * * * * * * * * * * 
subroutine ft:3( wold, px , alpha, xlarn ) 

C *** * * * X- ****** * * * * * * * * * * * * * *** * * * * * * * * * * * * * * * * * ****** * * * * * * * * * 

real *4 px < 25, 25 >, xlamt 25, 25 ), xold< 25, 25 ) 

camp ]. ex 2 0 , 2 1,22,2 3 , 4, 2 5 , z 6 , 7 

common n , n2 , n2 1 , n4 , n 1 , mn , mn2, nm3 , mn4 , m \ , m2 

x n 2 - n 4 

pin = 3 . 0 * a t an ( 1 . 0 ) / f 1 oat ( n ) 

xOO = 0,0 

>:22 = 0.0 

yOO = 0.0 

y22 = 0 . 0 

do 100 i = ml, m2 

do 100 j= ml, m2 

x 00 = x 00+px (i , j ) 

x22 = x 22-’-px ( i , j ) * ( ( — 1.0) ** ( i +j -mn4) ) 
yOO = y00+ xol d ( :i , j ) 

y22 = y 22+x ol d ( i , j ) * ( ( - 1 . 0 ) ** ( i h j -mn4 ) ) 

1 00 x 1 am ( i , j > = 0.0 

x 00 = 1 - 0 / y 00+ ( 1 - 0 - a 1 p h a ) * x 00 
xOO = l.O/xOO 

x 22=1. 0 / y 22 + ( 1 . G-a 1 ph a ) *x 22 
x 22 = 1,0/ X22 


do 150 m = ml, m2 
if ( m . 1 t . mn2 ) go to 985 
nnl = mn2 

go to 9B6 

935 nnl = mn3 

9S6 continue 

do 150 nn - nnl, m2 

150 1 am (m, nn ) = xOO + x22 * ( ( -1.0) **( m+nn-ir.n4) ) 

x 0 = p x ( m n 2 , m n 2 ) 

yO : = xol d (mn2, mn2) 

zO = cmpl x (cos (pi n ) , sin (pin)) 

2 1 ( 1 . 0 , 0 . 0 ) 

z 6 — c nip 1 ;; ( c 0 s ( p :i. n >< n 1 ) , —si n ( p i r 1 K n 1 ) ) 
z 7 = (1.0, 0.0) 




d o 

300 

k = 1, n 2 1 



= zl 



x 1 

- 0. 

0 


yi 

— o . 

0 


do 

200 

rn = mn3 ? m2 


X X 

== real (z2) 


xl 

= xl 

+ px ( rn, mn2 ) *xx 


yi 

= yi 

•h xold ( m, mn2 ) * 

200 

“T O 



* z 1 


Z o 

= ( 

1 - 0 , 0 . 0 ) 


do 700 1 = 1 , n 

it C k*l . eq. 1 ) go to 700 

it ( k m eq. 1 • and. 1 - gt . n21 ) go to SOI 

i t ( k. eq. n21 . and. 1 . gt . n2 ) go to S00 

x 2 = 0-0 
y 2 = 0.0 
z4 = z 7 

do 400 m ml, m2 
sb = ::3 

do 300 nn = mn3, m2 
x x = real ( £ 4* ) 

x 2 = >: 2 + px ( m , nn ) * x x 

y2 = y2 +• xold( m 9 nn ) k 
300 25=25*23 

it ( rn . e q . m2 ) g o t o 4 0 0 
z4= z4 ^ zl 
400 cant i nue 
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x k 1 = x 0 +■ 2 . 0 * ( x 1 + x 2 ) 
ykl - yO + 2.0 * ( yl + y2 ) 
x k 1 = 1.0 / y k 1 + ( 1.0 ~ alpha ) * >; k 1 

x k 1 = 1.0 / xkl 

2 2 = z 7 

do 600 rn = m3 , m2! 

it ( m.lt.mn2 ) go to 987 

z 5 = ( 1 . 0 y 0.0 ) 

nn 1 - mn2 

go to 98B 

987 z5 = z 3 
nn 1 := m n 3 

988 continue - 

do 500 nn = nnl , m2 
it (k*l.eq.n21 ) go to 989 

xlam (ir>,nn> = >: 1 am ( rn H nn ) +2 . 0*x k 1 * real (z2*z5) 
go to 500 

989 x 1 am (m,nn! =x 1 am ( rn ,nn) h-x k 1 *r ea 1 ( z 2*:: 5 ) 


500 

_ c 

aZ i_J 

z 5 * 

Z o 

□ 00 

z 2 - 

2 2 * 

z 1 

700 

M — 

z 3 * 

z 0 

301 

1! 

•H 

N 

z 1 * 

zO 


200 


300 


400 


994 


2 0 

- c mp 1 x ( c os (pin), sin (pin) ) 

2 1 

-• ( i . 0 , 0 .. 0 ) 


2 6 

= cm p 1 x < cos (pin *n 1 ) , -si 

n (pi 

2 7 

- ( 1 . 0 , 0 . 0 ) 


do 

800 k ~ 1, n21 


-r '“> 

- zl 


:•! 1 

= 0.0 


do 

200 m - mn3, m2 


x 1 

= xl + x 1 am ( m, mn2 > *real ( 

22) 


- 22 * zl 


^ 3 

= ( 1 . 0 , 0 . 0 ) 


do 

700 1 = 1 , n 


if 

( kxl.eq.l ) go to 700 


if 

(k . eq. 1 . and . 1 . gt. . n21 ) go 

to 

if 

( k . eq. n2 1 . and . 1 . gt . n2 ) 

go t 

X 

= 0 . 0 


2 4 

= z7 


do 

400 rn = ml, m2 


2 5 



do 

300 nn = inn 3, m2 


X 2 

= x2 + x 1 am (m, nn ) * real ( 

2 4* 

2 5” 

: z5*z3; 


if 

( m. eq. m2 ) go to 400 


2 4 

= z 4 * z 1 


cont i nue 


xkl 

= x 0 + 2.0 # ( x 1 +x 2 ) 


if 

C xkl.gt.0.0 ) go to 994 


_L . . _ r-v i—\ r—, 

format ( 5x , bad solution, 

f ( 


j 01 


*( CHECK POINT NO- 2 ( SUBROUTINE FT)) 

stop 

conti n ue 

= a m i n 1 ( x rn i n 1 , >i k I ) 
k m i n 2 — a max 1 ( xmin2 P xkl ) 
x k 1 = 1. 0 / x k I 


0 

) 


22 = 27 

do 600 m~ ml, m2 

.if ( rn.lt.mn2 ) go to 995 

z 5 = ( 1.0, 0.0 ) 

nn 1 = mn2 

go to 996 

995 z 5 = 23 
nn 1 = mn3 

996 continue? 

do 5 0 0 n n = n n 1 , m2 
if* ( k * 1 . 0 q . n 2 1 ) go to 997 

p x ( m , nn ) = px ( m, nn ) +2.0 * xkl * real ( 

cjo to 500 

p v ( m „ n n ) = px ( m , nn ) + 2-^s j ) 
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22*25 ) 


997 
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500 

Zu = 2 5 

* z3 

600 

7 o — 

* zl 

700 

v r: 2! 3 

* z0 

BO 1 

Zl = Zl 

* z 0 


z7 = z6 

* con j g ( z 4 ) 

BOO 

conti nue 


rJa 900 rn= ml , m2 
d o 900 n n = m 1 , rnn 2 

900 

p x ( m j, n n 

) -■ p ( rn 1 + m 2 - m , m 1 -h m 2 


px0~ px ( 

mn2, mn2 ) 


x x = p x 0 

/ xn2 


do 1000 

m = ml y m2 


do 1000 

nn = ml, m2 


px ( m, nn 

) = px ( m ? nn)/px0 

1000 

>: 1 am (m, nn ) = 1 am ( m, nn ) *>: >: 


x m i n 2 = 

1.0 / ( mi n2 *■ >:>: ) 


x mini = : 

: /it i n 1 * >; 


return 

end 


?-■.* 


;-.VirtV ■~**i~-**j***» mtm*. r* —■*■». -*• > 
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